Many times, linearity between a dependent variable and a set of predictors is not a suited assumption. Normal linear models rely on a normality assumption, and this may be a limitation when the dependent variable is discrete, e.g. with support \(\{0,1\}\) or \(\{0,1,2, \ldots, \}\), or it is continuous and its range spans from 0 to \(+\infty\). In normal linear models, the main assumption is (including the intercept \(x_{ij} = 1\)): \[ E[Y_i|{\bf x}_i]= \eta_i = \sum_{j=1}^{p}\beta_jx_{ij}.\] Instead, the generalized linear model consider modelling a suitable transformation of the conditional expectation by using covariates for a broad range of distributions, that is (including the intercept \(x_{ij} = 1\))
\[g(E[Y_i|{\bf x}_i])= \eta_i =\sum_{j=1}^{p}\beta_jx_{ij}.\] The main ingredients of a GLM are: the distibution of \(Y_i\), which belongs to the exponential dispersion family, the linear predictor (\(\eta_i\)), and the link function \(g(\cdot)\) that must be twice differentiable.

The following example is from Chapter 5 of Data analysis using regression and multilevel/hierarchical models, Gelman, A., & Hill, J. (2007), Cambridge University Press.

It is known that water in some wells in Bangladesh and other South Asian countries might be contaminated with natural arsenic. The effect of the arsenic ingestion on health is a cumulative poisoning, and the risk of cancer and other diseases is estimated to be proportional to the dose intake.

A research team from the United States and Bangladesh examined the water of several wells in Araihazar district of Bangladesh and they measured the arsenic level classifying the wells as safe if the arsenic level was below the 0.5 in units of hundreds of micrograms per liter, or unsafe if it was above 0.5. All the people using unsafe wells have been encouraged to switch to a nearby well. Indeed the contamination does not depend on the proximity and close wells can have very different levels of arsenic.

The dataset in the file Wells.csv contains information on whether or not 3020 households changed the unsafety wells after few years from the arsenic investigation. The variables are

• switch: whether or not the household switched to another well from an unsafe well: no (0) or yes (1);
• arsenic: the level of arsenic contamination in the household’s original well (in hundreds of micrograms per liter); note that all are above 0.5, which was the level identified as “unsafe”;
• distance: distance to the closest known safe well (in meters);
• education: education of the head of the household (in years);
• association: whether or not any members of the household participated in any community organizations: no (0) or yes (1).

Goal: Using a logistic regression, we would determine which are the main factors of well switching among the users of unsafe wells.

We start by reading the data, by analysing the structure and exploring some characteristics of the data.

We start by modelling the effect of the variable distance on the variable switch. Since the outcome variable is a binary variable, the model chosen is the logistic regression model. \[ \rm{Pr}(Y_i=1| dist_i)=E(Y_i | dist_i)=p_i \quad \quad \text{logit}(p_i)= \eta_i = \beta_1 +\beta_2 dist_i\] where \(\text{logit}(x)=\log(x/(1-x))\) is a function that maps the range \((0,1)\) to the range \((-\infty,\infty)\).

The intercept (\(\beta_1\)) can only be interpreted assuming zero values for the other predictors. Here, when dist\(=0\) the estimated probability of switching for a family who lives close to a safe well is about \(65\%\), given by \[\hat P(Y_i = 1 | dist_i = 0) = \frac{e^{\hat \beta_1}}{1+e^{\hat \beta_1}}= \frac{1}{1+e^{-\hat \beta_1}}\]

Interpreting the coefficients for the predictor dist on the probability scale: we can evaluate the difference of the logistic regression function by adding 1 to the predictor dist. This difference corresponds to a difference on the probability scale but requires the choice of the specific value of the predictor, that is \[P(Y_i = 1 | dist_i = c+1) - P(Y_i = 1 | dist_i = c) = \frac{1}{1+\exp(-(\beta_1 + (c+1) \times\beta_2))} - \frac{1}{1+\exp(-(\beta_1 + c \times\beta_2))}\] The mean of the predictor is a good starting point since it corresponds to the steepest point of the inverse logit function. Thus, adding 1 to dist (that is, adding 1 meters to the distance to the nearest safe well) corresponds to a negative difference of about \(0.15\%\) in the estimated probability of switching.

Of course, you can evaluate the (estimated) probability difference between two different distances.

Note that the coefficients here does not have a linear effect on the probability that the outcome is 1 because of the non linearity of the model. The curve of the logistic function requires us to choose where to evaluate changes, if we want to interpret on the probability scale. Below the difference on the (estimated) probability of switching well due to an increase of one unit on the predictor from the 99th percentile

The effect of one unit increase of the variable dist on the inverse logit seems low, but this is misleading. Indeed, the distance is measured in meters, so this coefficient corresponds to the difference between, say, a house that is 48 meters away from the nearest safe well and a house that is 49 meters away. To improve interpretability of the results we rescale distance in 100-meter units. Thus, adding 100 meters to the distance to the nearest safe well corresponds to a negative difference in the probability of switching of about \(15\%\).

• Note that \(\hat \beta_2\) of the new model is 100 times the \(\hat \beta_2\) of the original one

• Also note the following equivalence: the difference on probability (of switching well) scale for a unit increase of the new model corresponds to 100 times the same quantity in the original model

We consider adding the arsenic level in our model. We expect that the higher the arsenic concentration the higher the probability of switching well, that is, translated in model terms, a positive sign for the arsenic coefficient. We plot the histogram of arsenic concentration before fitting the model and analysing the results.

Basically, we fitted the model \[\text{logit}(P(Y_i|dist^{*}_i, arsenic_i))= \eta_i = \beta_1 +\beta_2 \texttt{dist}^{*}_i + \beta_3 \texttt{arsenic}_i\]

(note \(dist^* = dist/100\))

• Given two wells with the same arsenic level (notice that it cannot be 0!), the estimated logit probability of switching decreases of \(.9\) every 100 meters in distance to the nearest safe well. \[\text{logit}(P(Y_i|dist^*_i=c+1, arsenic_i = k)) - \text{logit}(P(Y_i|dist^*_i=c+1, arsenic_i = k)) = \beta_2\]
• For two equally distant nearest wells, a difference of 1 in arsenic concentration corresponds to a 0.46 positive difference in the logit probability of switching. \[\text{logit}(P(Y_i|dist^*_i=c, arsenic_i = k+1)) - \text{logit}(P(Y_i|dist^*_i=c, arsenic_i = k)) = \beta_3\]

The following plots show the fitted model with two predictors. Numerical interpretation of one predictor effect on probability scale needs to be done choosing a value for the other predictor.

• Left: probability of switching well vs distance when arsenic is equal to 0.5 and 1.5

• Right: probability of switching well vs arsenic when the distance is equal to 0 and 50

The following plots show the fitted model with two predictors. Numerical interpretation of one predictor effect on probability scale needs to be done choosing a value for the other predictor.

• Left: probability of switching well vs distance when arsenic is equal to 0.5 and 1.5

• Right: probability of switching well vs arsenic when the distance is equal to 0 and 50

In general terms, the fitted model can be summarized as: switching is easier if there is a nearby safe well, and if a household’s existing well has a high arsenic level, there should be more motivation to switch.

We now want to verify if there is a statistically significant effect of the interaction between the distance and the arsenic concentration on the probability of switching well.

We fix the predictors to their mean and so the estimated probability of switching well for a averaged distance and averaged arsenic concentration is about 0.59. \[\hat P(Y_i = 1|dist^*_i = \overline{dist^*}, arsenic_i = \overline{arsenic}) = \frac{\exp(\hat \beta_1 + \hat \beta_2 \overline{dist^*} + \hat \beta_3 \overline{arsenic} + \hat \beta_4 \overline{dist^*} \times \overline{arsenic})}{1+\exp(\hat \beta_1 + \hat \beta_2 \overline{dist^*} + \hat \beta_3 \overline{arsenic} + \hat \beta_4 \overline{dist^*} \times \overline{arsenic})}\]

At the mean level of arsenic in the data, each 100 meters of distance corresponds to an approximate 22% negative difference in probability of switching. \[\hat P(Y_i = 1|dist^*_i = c + 1, arsenic_i = \overline{arsenic}) - \hat P(Y_i = 1|dist^*_i = c , arsenic_i = \overline{arsenic})= \] \[=\frac{\exp(\hat \beta_1 + \hat \beta_2 (c + 1) + \hat \beta_3 \overline{arsenic} + \hat \beta_4 (c+1) \times \overline{arsenic})}{1+\exp(\hat \beta_1 + \hat \beta_2 (c+1) + \hat \beta_3 \overline{arsenic} + \hat \beta_4 (c+1) \times \overline{arsenic})}- \frac{\exp(\hat \beta_1 + \hat \beta_2 c + \hat \beta_3 \overline{arsenic} + \hat \beta_4 c \times \overline{arsenic})}{1+\exp(\hat \beta_1 + \hat \beta_2 c + \hat \beta_3 \overline{arsenic} + \hat \beta_4 c \times \overline{arsenic})}\]

At the mean level of distance in the data, each additional unit of arsenic corresponds to an approximate 11% positive difference in probability of switching (the mathematical expression is similar to the expression above)

The interaction term can be seen (on the logit scale): as the effect (-0.18) added to the coefficient for distance for each additional unit of arsenic; indeed

\[{\rm logit}(P(Y_i|dist^*_i, arsenic_i = k)) = \beta_1 + \beta_2 dist^*_i + \beta_3 k + \beta_4 dist^*_i \times k\]

\[{\rm logit}(P(Y_i|dist^*_i, arsenic_i = k + 1)) = \beta_1 + \beta_2 dist^*_i + \beta_3 (k + 1) + \beta_4 dist^*_i (k + 1) = \] \[ = \beta_1 + (\beta_2 + \beta_4) dist^*_i + \beta_3 (k + 1) + \beta_4 dist^*_i \times k \] So, the importance of distance as a predictor increases for households with higher existing arsenic levels (\(\hat \beta_2 + \hat \beta_4 = -0.58 - 0.18\))

Further, \[ {\rm logit}(P(Y_i|dist^*_i, arsenic_i = k + 1)) - {\rm logit}(P(Y_i|dist^*_i, arsenic_i = k)) = \beta_3 + \beta_4 dist^{*}_i\] So, at the average distance the difference in estimated logit probabilities is 0.46

The interaction term can be seen (on the logit scale): as the effect (-0.18) added to the coefficient for arsenic for each additional 100 meters of distance; indeed

\[{\rm logit}(P(Y_i|dist^*_i=c, arsenic_i)) = \beta_1 + \beta_2 c + \beta_3 arsenic_i + \beta_4 c \times arsenic_i\]

\[{\rm logit}(P(Y_i|dist^* = c+1, arsenic)) = \beta_1 + \beta_2 (c+1) + \beta_3 arsenic_i + \beta_4 (c+1) arsenic_i = \] \[ = \beta_1 + \beta_2(c+1) + (\beta_3 + \beta_4) arsenic_i + \beta_4 c \times arsenic_i \] So, the importance of arsenic as a predictor decreases for households that are farther from the existing safe wells (\(\hat \beta_2 + \hat \beta_4 = 0.56 - 0.18\))

Further, \[ {\rm logit}(P(Y_i|dist^*_i=c+1, arsenic_i = k)) - {\rm logit}(P(Y_i|dist^*_i=c, arsenic_i = k)) = \beta_2 + \beta_4 arsenic_i\] So, at the average value of arsenic the difference in estimated logit probabilities is -0.88.

The following plots show that the differences in switching associated with differences in arsenic level are larger if you are close to a safe well, but with a diminishing effect if you are far from any safe well. Comparing with the same plot without interaction, note that the two curves mainly differ for higher values of distance.

The following plots show that the differences in switching associated with differences in arsenic level are larger if you are close to a safe well, but with a diminishing effect if you are far from any safe well. Comparing with the same plot without interaction, note that the two curves mainly differ for higher values of distance.

The model selection can be carried out also by using the AIC. It judges a model by how close its fitted values tend to be to the true expected values, as summarized by a certain expected distance between the two. Because smaller AIC is better, the AIC penalizes a model for having many parameters.

Outcome variable in the logistic regression can take only 0 and 1 values. The plot of the residuals defined as \[{\rm residual}_i = y_i - \hat {\pi}_i = y_i-\text{logit}^{-1}({\mathbf x}^\top_i\hat \beta)\] versus the fitted values (\(\hat \pi_i= \text{logit}^{-1}({\mathbf x}^\top_i\hat \beta))\) is not useful, that is considering \((\hat \pi_i, y_i - \hat \pi_i)\), because the points belong to two parallel lines with slopes -1. Let’s consider the model labelled as fit4.

Outcome variable in the logistic regression can take only 0 and 1 values. The plot of the residuals defined as \[{\rm residual}_i = y_i - \hat {\pi}_i = y_i-\text{logit}^{-1}({\mathbf x}^\top_i\hat \beta)\] versus the fitted values (\(\hat \pi_i= \text{logit}^{-1}({\mathbf x}^\top_i\hat \beta))\) is not useful, that is considering \((\hat \pi_i, y_i - \hat \pi_i)\), because the points belong to two parallel lines with slopes -1. Let’s consider the model labelled as fit4.

The following plot shows the binned residuals vs the the estimated probability of switching. Grey lines depict \(\pm 2\) standard error (se) bounds, so we expect that 95% of the binned residuals falls between these two lines. The se is defined as \(\sqrt{p(1-p)/n}\), where \(n\) is the number of data used to compute each average residual.

The following plot shows the binned residuals vs the the estimated probability of switching. Grey lines depict \(\pm 2\) standard error (se) bounds, so we expect that 95% of the binned residuals falls between these two lines. The se is defined as \(\sqrt{p(1-p)/n}\), where \(n\) is the number of data used to compute each average residual.

A pattern can be identified looking at the plot of the binned residuals with respect to the arsenic level. Indeed, there are large negative binned residuals for some households with wells with small values of arsenic concentration. These points correspond to people that are less likely to switch well with respect to what predicted by the model (almost \(-20\%\)). Moreover, the distribution of the residuals shows a slight pattern: the model seems to underestimate the probabilities of switching for medium arsenic level values while overestimates for high arsenic concentrations. So we plot the binned residuals vs the the predictors.

Error rate is defined as the proportion of cases for which the fitted probabilities are above a threshold (commonly 0.5) and the outcome value is 0 or the fitted probabilities are below the specified threshold and the outcome value is 1, that is when the prediction-guessing is wrong: \[er = \frac{1}{n}\sum_{i=1}^{n}((y_i=0 \, \& \, \hat p_i >0.5)\, | \, (y_i=1 \, \& \, \hat p_i <0.5))\] It should be a value between 0 (perfect prediction-guessing) and 0.5 (random prediction-guessing).

Usually we expect the error rate to be lower than the error rate of the null model, that is the model including only the intercept. The estimated proportion for the null model will be \(\hat p=\sum_i y_i/n\), that is the proportion of 1’s in the data. The error rate for the null model is \({\rm min}(p, 1-p)\), and it corresponds to assign 1 (or 0) to all the fitted values.

So if our model improve the prediction-guessing of a simple logistic regression, we expect the error rate of our model to be lower than the error rate for the null model. Our last fitted model (model labelled as fit4), with an error rate of \(0.38\), correctly predicts \(62\%\) of the switching choices, determining only a \(4\%\) of improvement if compared to the simply guessing that all the households will switch.

Consider the error rate as the equivalent of the \(R^2\) for the linear model, it can be a useful measure of the model fit but can’t substitute a deeper evaluation of the model looking at coefficients significance and at the diagnostic plot. Analysing error rate means to extract information from the confusion matrix. Note that we can vary the threshold 0.5, and better assessment can be done by analysing the ROC curve.